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I. INTRODUCTION 


The possibility that a large asteroid may impact the Earth delivering an energy in 
excess of 30 MT of TNT is a very real threat. One only has to look at the Moon in the 
night sky to gain an appreciation for the magnitude and probability of such an impact. 
How the human race deals with the threat of such an impact is of significant interest to the 
planet Earth as a whole. 

As discussed by Rather et. al. (1992), as early as 1705, when Edmund Halley 
wrote A Synopsis of the Astronomy of Comets, there has been speculation that an 
extraterrestrial object might impact the Earth. the possibility of such impacts was not 
perceived as a threat to life on Earth until the late 1940's when the Moon's craters were 
fully understood to be the result of impact events, not volcanism, and that the Earth is 
subject to the same impact hazard. Rather et. al (1992) also mention that the magnitude 
of the threat was not fully realized until the late 1970s to early 1980s. The perceived 
threat began to stand upon a solid foundation with the publication by Alvarez et. al. 
(1980) ofa theory on the extinction of the dinosaurs due to the impact of a large asteroid 
or comet 65 million years ago. The theory put forth by Alverez et. al. suggests that an 
impact by a 10 km diameter asteroid at the Chicxulub site off of the Yucatan peninsula 
indirectly caused the extinction of 60% of all life on the Earth, including the dinosaurs. 

Since the awareness of the possibility of an asteroid impacting the Earth began 
developing in the mid 20" century, there have been several "near misses" recorded. 
Perhaps the most spectacular near miss was a large fireball created by an object racing 
across the daytime sky in a northerly direction that entered the Earth's atmosphere above 
Utah in the United States and exited above Alberta, Canada on the 10" of August, 1972. 
The object observed was determined to be an asteroid upwards of 30 meters in diameter 
as reported in Sky and Telescope Magazine (1972). Had this asteroid's trajectory been 
ever so slightly different, mankind would have had its first opportunity to observe a large 
object impact the Earth. Such an impacting object would carve out a crater 200 to 500 
meters in diameter. Since then, several other asteroids and comets have been detected 
passing by the Earth at distances less than a few hundred thousand kilometers. Apollo 
Asteroid 1989FC, referenced by the AIAA Space Systems Technical Committee; Asteroid 


1991BA, reference by Scotti et. al. (1991); and Asteroid 1996JA1 referenced by Jaroff 
(1996) are three asteroids discovered recently passing very close to Earth. Onan 
astronomical scale, these close approaches are essentially impacts. 

In April of 1990 the American Institute of Aeronautics and Astronautics issued the 
position paper entitled Dealing with the Threat of an Asteroid Striking the Earth that 
briefly described the implications of an asteroid impact. The AIAA found that "Earth- 
orbit-crossing asteroids clearly present a danger to the Earth and its inhabitants." It was 
recommended that a "systematic and open program" for detection of Earth crossing 
asteroids be established as well as a study to "define systems which can deflect or destroy, 
or significantly alter the orbits of, asteroids predicted to impact the Earth." A few search 
programs existed prior to the position paper and a few others have begun subsequently. 

As a result of the awareness of the possibility of an asteroid or comet impacting 
the Earth, the National Aeronautics and Space Administration has sponsored two 
workshops to study the fundamentals of the impact and impact mitigation problem. The 
first workshop is summarized in the report entitled "The Spaceguard Survey: Report of 
the NASA-International Near-Earth-Object Detection Workshop." The second workshop 
is summarized in the report entitled "Near-Earth-Object Interception Workshop." The 
concentration of the workshops have been related to assessing the magnitude of the threat, 
impact effects and hazards to the Earth, as well as the political implications of developing 
an impact mitigation capability. Several books have been published on the matter, both 
technical and non-technical, and more than one Hollywood movie has been based on the 
subject. Two spacecraft exploration missions, Near Earth Asteroid Rendevouz (NEAR) 
and Clementine, have included intercepts of asteroids as a major part of their mission to 
study the nature of asteroids. However, little "astrodynamical analysis" has been 
performed on the feasibility of the impact mitigation problem. 

The astrodynamical analysis of feasibility is where this thesis is intended to fit into 
the larger problem. Presented is an analysis of the impact and mitigation problem based 
on a two-dimensional and two-body analysis. It is intended to be a first order 
approximation for the solution of a larger problem. However, it is a more rigorous 


treatment of the astrodynamics of the hazard than previous published analyses, such as 


Ahrens and Harris (1994). Included in the following analysis is the periodic nature of the 
problem as well as the near term effects. The analysis is not concerned with object 
detection or orbit prediction, but instead centers on how impulses applied to an asteroid at 
various points on the asteroid's orbit affect the outcome when there is a presumption of 
collision otherwise. Mission design for mitigation is not a goal of the following analysis; 
however, the analysis tool presented may be utilized in determining a first order estimate 
for optimizing the time and position of asteroid intercept for impact mitigation. 
Presented first is an astronomical development of the existence of an asteroid 
impact problem from the origins of the solar system to the record of past impacts on the 
Earth followed by a brief description of two-body orbits. The problem is then presented 
with governing assumptions and a method of solution. The solution method is then 
assessed and validated followed by an analysis of an asteroid impact scenario with a 
discussion of the results. The analysis method is then applied to the asteroid Toutatis 


which will make several close approaches to the Earth over the next decade. 
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Il. SOLAR SYSTEM OBJECT ORIGINS 


A. SOLAR SYSTEM FORMATION 

Current theories of formation of the Solar System stem from the post initial 
expansion universe environment of gas, dust, radiation, and magnetic fields in a 
nonuniform distribution. Bouyed by the interplay of gravitational, magnetic, and pressure 
forces local mass concentrations began to form. This interplay of forces on the non- 
uniform environment sets up initial angular momentum conditions for very large three- 
dimensional structures. These large rotating structures began to coalesce into accretion 
disks around central masses. These central masses eventually reached critical mass for 
fusion to occur and stars emerged. We call our local star the Sun. 
B. PLANETARY FORMATION 

In a similar fashion to the stellar formation, the accretion disk around the Sun 
provided the environment for smaller mass concentrations and accretion disks to form thus 
producing small scale structures called planetismals composed of the basic chemical 
elements in varying quantities. Gravitational attraction and relative motion of these 
planetismals caused them to collide with one another and cohesive forces enabled some to 
remain attached forming larger structures. After enough of these interactions occurred, 
the planets began to form. Unlike for the stellar conditions, the planets do not possess the 
critical mass to initiate and sustain a fusion reaction. This allows for large scale assembly 
of solid, aqueous and gaseous structures to take place in quantities and composition 
proportional the relative percentages of chemical elements present in the planetismals. 
These structures are more familiar on the Earth as the crust, the oceans and the 
atmosphere. Analysis of the elements present from the impact delivery mechanism and 
current known compositions of asteroids and comets suggests that this was the mechanism 
of organic and non-organic material delivery that formed the Earth as proposed by Chyba 
et. al. (1994). 
Cc, ASTEROID AND COMET FORMATION 

However, not all of the accretion disk surrounding the Sun coalesced into either 
the Sun, the Planets or their satellites. The remainder of the planetismals continue to be 


dispersed throughout the Solar System in the form of asteroids and comets. The asteroids 


being located primarily within the inner Solar System, and the comets existing in the outer 
Solar System and beyond. 

The asteroids are mainly concentrated in the asteroid belt located between the 
orbits of Mars and Jupiter. The remainder of the asteroids are dispersed throughout the 
solar system in elliptical orbits lying mainly inside the orbit of Jupiter. There are several 
theories of how the asteroid belt came to exist. These ideas include the destruction of a 
planet, or the inability of a planet to form, due to the tidal forces of Jupiter. There are 
several theories accounting for the formation of the asteroids not located in the asteroid 
belt. These ideas range from planetismals that have never collided and attached 
themselves to another planetary body, to cast off remnants of massive collisions of 
planetary bodies with very large planetismals, to objects from the asteroid belt that may 
have been perturbed by a passing object into a smaller orbit. All of the ideas have some 
amount of merit that give them validity. 

Comets are believed to originate from the Oort Cloud of cometary material 
orbiting the Sun at a distance of some 50,000AU. From this cloud, comets are believed to 
be injected into the solar system by orbital perturbations due to the gravitational field of 
passing stars. Once injected into the solar system, gravitational encounters with the Sun 
and Jupiter may further perturb the comets’ orbit and either “capture” the comet, so it 
remains within the solar system, or “assist” the comet on its way to interstellar space. 
Additionally, there is believed to exist a band of icy objects that extends from the orbit of 
Neptune at 30 AU out to as much as 100 AU. These objects are said to be located in the 
Kuiper Belt, so named after Gerard Kuiper, who first proposed their existence in 1951. 
The objects that form the Kuiper belt range in size and orbital characteristics to the extent 
that it is believed that Pluto may actually be a member of this group as discussed by Jewitt 
and Luu (1995). 

D. ASTEROIDS 


1. Location 
The asteroids are of particular interest as potential impact hazards in that they have 
a greater mass density than comets and are more likely to reach the surface of the Earth in 


a given impact scenario. The majority of the asteroids are located in the asteroid belt 


between the orbits of Mars and Jupiter. Very few asteroids have been detected beyond the 
orbit of Jupiter. This lack of detection may only be an effect of the limited capability of 
current detection sensors. However, a significant number of asteroids in orbits smaller 
than that of the typical asteroid belt object have been identified. These are the objects of 


primary concern for the problem of mitigation. 


Z Quantities 

Literally thousands of asteroids have been observed and identified orbiting the Sun. 
Of those, more than 300 are considered near Earth asteroids (NEA’s) and pose a threat as 
a potential impacting object. Of greater concern are the subset of NEA’s dubbed Earth 
crossing asteroids (ECA’s) which are currently about 200 in number. The orbits of the 
ECA’s are such that they allow for the possibility of impact with the Earth at some future 
date. These ECA’s range in size from 10 km down to 0.1 km, which is currently the limit 
of detection of asteroids by Earth based sensors. Estimates for ECA’s near and above 10 
km in diameter indicate that all of the asteroids have been identified and that for the ECA’s 
of 1 km in diameter or less, the identified asteroids represent only about 10% of those that 
are believed to exist, as stated by Grieve and Shoemaker (1994). 


oe Classification 

Classification of the ECA's are determined with respect to the Earth's orbital 
extrema. The Earth's perihelion and aphelion distances are 0.9833 and 1.0167 AU 
respectively. The orbits of the ECA’s all have perthelia less than the aphelion of the Earth 
orbit and aphelia greater than the perihelion of the Earth orbit. The ECA’s have been 
subdivided into three classes, the Atens, Apollos and Amors, based on their orbital 
characteristics as discussed by Rabinowitz et. al. (1994). Table 1 summarizes the 
distinction between the three classes. Figures 1, 2 and 3 depict typical orbits of each of 


the three classes. 
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4, Physical Properties 

The physical properties of the asteroids are generally that of “rocky”, regularly 
shaped spinning objects. Estimates of asteroid densities range from a more cometary 
density of 2x10° kg m” to a dense asteroid of 5x10° kg m” with a mean density of about 
3x10° kg m’. Asteroids have been observed that are composed of a single solid mass as 
well as multiple mass centers either physically connected or gravitationally bound at a 
contact surface. As stated by Winters (1996), some of the asteroids may actually be 
aggregates of numerous smaller bodies that are gravitationally bound together. This 
hypotheses is further supported by analysis of object spin motion. It appears the many 
asteroids spin at angular rates sufficiently slow as to permit gravitational binding. The 
case of the comet Shoemaker-Levy 9 supports these theories in that it was gravitationally 
separated by tidal forces from a previous passage of Jupiter prior to its 1994 impact. 
E. COMETS 


1. ~° Location 

Early in the 20" century, Jan Oort hypothesized that virtually all of the comets 
originate in a cloud of cometary matter beginning some 50,000AU from the Sun. The so 
called Oort cloud is assumed to be a uniformly distributed spherical shell that surrounds 
the solar system and extends to approximately half of the distance to the nearest star, 
Alpha Centauri. From this cloud, comets are injected into the solar system by orbital 
perturbations of passing stars. Once inside the solar system the comets may be further 
perturbed by the planets. The planetary perturbations may “capture” the comet in the 
interior of the solar system or it may assist the comet on its way ejecting it from the solar 


system for all time. 


Ze Quantities 
The number of comets in the Oort cloud is believed to be diminishing, however the 
remaining number is believed to be on the order of 10'*. The quantity of comets that exist 


within the solar system is far smaller, only 100 or fewer have been identified. 


a Classification 
The “captured” comets are classified as periodic comets and are subdivided into 


two groups as summarized by Shoemaker et. al. (1994). Members of “Jupiter’s Family” 
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have aphelia close to Jupiter’s mean orbital distance from the Sun and have periods of less 
than 20 years. Members of the “Halley family” are the so called long period comets and 


have orbital periods greater than 20 years but less than 200 years. 


4. Physical Properties 

Comets are primarily composed of “rock” and “ice.” As such they have been 
sometimes called “dirty snowballs” or “icy dirtballs’” depending on their relative 
composition. The icy material is believed to be composed of water, methane or ammonia. 
The rocky material is a variety of carbonaceous substances. Their mean density is less 
than that of the asteroids and is estimated to be about 1000 to 2000 kg m”. As a comet 
“burns” off its icy material from repeated encounters with the Sun, the nucleus may in fact 


become an asteroid of a somewhat smaller dimension. 


1] 





Il. IMPACT RECORD 


Occasionally an asteroid or a comet’s orbit is such that it actually hits another 
object within solar system such as the event, widely celebrated in the media, of comet 
Shoemaker-Levy 9 impacting the planet Jupiter on the 16" of July, 1994. The very same 
mechanism that caused the solar system, in particular the Earth, to form and sustain life 
now threatens the very same life. 

A. LUNAR IMPACT RECORD 

One theory of the origin of the Moon, as mentioned by Chapman and Morrison 
(1994) describes it forming as the result of a Mars size object impacting the Earth early in 
its existence . The ejecta from the event is believed to have behaved in such a fashion as 
to remain in orbit and coalesce into what is now the Moon. Further evidence of such 
impacts is the cratered face of the lunar surface. The Moon displays literally millions of 
impact sites. Since the Moon has no atmosphere or large scale geologic processes there is 
no erosion of the impact record, unless by subsequent impact the previous impact 
structure is destroyed. Thus, the Moon serves as a reminder for all time of the nature of 
the impact hazard. 

B. TERRESTRIAL IMPACT RECORD 

On our planet Earth, the atmosphere, oceans, volcanism and plate tectonics tend to 
erode past impact sites. As discussed by Grieve and Shoemaker (1994) there are currently 
about 140 known impact craters around the world. The locations of these sites are 


displayed in Figure 4. 
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Figure 4. Known Impact Site Locations from Grieve and Shoemaker (1994) 


These impact craters range in size from 1.5 km in up to 200 km in diameter. 
Perhaps the best example of a classic impact crater is the 1.5 km diameter Barringer Crater 


located in Arizona, see Figure 5. Barringer Crater is believed to be the most recent impact 
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site on Earth having formed around 50,000 years ago. 
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However, the most recent impact event to have resulted in surface damage is 


believed to be the Tunguska event which took place over northern Siberia in 1908. The 
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Tunguska event is believed to be the result of a 60 m diameter comet or light asteroid that 
exploded 10 km up in the atmosphere releasing some 30 MT of energy leveling 2500 
square kilometers of forest. The Tunguska event is characteristic of a small scale impact 
occurrance. On the large scale end of the impact spectrum hes the K/T impact event so 
called by its occurrence in time at the boundary of the Cretaceous and Tertiary periods. It 
is believed that an enormous asteroid some 10 kilometers in diameter created the 200 km 
diameter Chicxulub impact site beneath the Gulf of Mexico off the coast of the Yucutan 
peninsula and is responsible for the extinction of 60% of the living species on Earth 65 
million years ago. 

There have been numerous recent close calls of asteroids impacting the Earth. 
Jaroff (1996) describes a near impact as recently as June of 1996 where an object roughly 
600m in diameter passed within 450,000 kilometers (about 70 Earth radii) of the Earth. 
c. OTHER PLANETARY IMPACTS 

Numerous impact sites have been observed by planetary exploration spacecraft that 
have been sent to Venus and Mars. As mentioned above, Comet Shoemaker-Levy 9 
impacted the planet Jupiter in July of 1994. 
D. IMPACT SCALE 

The energies released by the impact of asteroids on the Earth are quite large. It is 
conventional to express the impact energy in terms of megatons of TNT (MT), where 
IMT = 4.2x10'°J. With the assumption of a mean asteroid density of 3000 kg m” a size 
versus impact velocity relation may be made for a given impact energy. Figure 6 displays 
an estimate of asteroid mass versus impact velocity for a range of impact energies from 
tens of megatons of TNT up to a petaton (10'°) of INT. The diameter estimate assumes 
an effective spherical radius corresponding to the mass for the same impact energy and 


velocity. 
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Figure 6. Impact Scale 


Estimates of impact energy and terrestrial devastation have lead to the 
classification of impacts as local, regional, and global events as described in Morrison et. 
al. (1994). The distinction between these is somewhat blurred, however it is proposed 
that local events correspond to impact energies in the vicinity of 30 MT or less which will 
affect approximately 0.001% of the Earth’s surface area or about the size of a large 
metropolitan area. Regional impact events are those that release energy in the vicinity of 
300 MT to 3x10* MT and affect about 0.1% of the Earth’s surface or about the size of a 
large state. Global events are considered to be impacts that release energies near and 
above 3x10° MT affecting approximately 10% of the Earth’s surface area or about the size 
of a large country. These large events may cause such disruption of the ecological 
environment by particulate injection into the atmosphere that greater that 25% of the 


Earth’s population may be eliminated. Figure 6 also shows several representative impact 
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scenarios for a typical impact velocity of 20 km s’. The Tunguska event is depicted by the 
'*' symbol. Tunguska represents a nearly classical local impact event of an object some 60 
m in diameter. The recent May 1996 near miss is depicted by the '+' symbol and 
represents a regional event for an object some 600 m in diameter. A global impact is 
indicated by the ‘o' symbol and corresponds to an asteroid some 1.5 km in diameter. The 
'®' symbol represents the K/T event that created the Chicxulub impact site and 
corresponds to a 10 km diameter object. It is noted that the K/T event is well above the 
threshold for global catastrophe. 


1% 
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IV. ORBITS 


A. CONIC SECTIONS 

In a simple inverse square gravitational field, orbits take the shape of conic 
sections. For an object that lacks sufficient energy to escape the gravitational attraction of 
the central body, the resultant orbit will take the shape of an ellipse. For an object with 
sufficient energy to escape the gravitational attraction of the central body, the resulting 
orbit will take the shape of an hyperbola. The transition from an elliptical to hyperbolic 
orbit occurs along an “escape” trajectory shaped as a parabola. Straight line orbits that 
lead to collision of the orbiting body with its mass center are special cases of elliptical, 


parabolic, and hyperbolic orbits. 


B. COORDINATE FRAMES 


1. | Three-dimensional 

To define the physical problem in time and space a reference frame needs to be 
established. For the general three-dimensional case a Cartesian Sun-centered inertial, or 
Heliocentric, coordinate frame is chosen. Figure 7 displays the orientation of the 


Heliocentric coordinate frame. 








first day 
first day af spring 


of summer 


first day 
of winter 
first day 
af autumn 
Xe (Seasons are for Northern Hemisphere ) 
vernal equinox 
direction 


bad fi 


Figure 7. Heliocentric Coordinate System from Bate et. al. (1971) 
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The direction in the ecliptic plane from the Sun to the First point of Aries serves as 
the primary coordinate direction in the Heliocentric frame. The ecliptic normal serves as 
the “vertical” coordinate and is defined as positive in the northern half of the celestial 
sphere. The last coordinate direction is defined by taking the cross product of the 


previous two coordinate directions. 


De Two-dimensional 

To simplify matters for the present analysis, a two-dimensional planar perifocal 
coordinate frame is chosen. The principal axis defining a two-dimensional elliptical orbit is 
the direction toward periapsis from the primary focus. It is from this primary axis that the 
true anomaly is measured in a counterclockwise sense. The secondary axis is normal to 
the primary axis in a right handed sense. The resultant coordinate system is shown in 


Figure 8. 
Q 


Ig , Perihelion 


Figure 8. Pertfocal Coordinate System 
C. ORBITAL ELEMENTS 


ie Three-dimensional] 

In the general three-dimensional case six orbital elements are required to define a 
particular location in space of an object in orbit. Those orbital elements are the semi- 
major axis (a), eccentricity (e), inclination to the ecliptic (i), longitude of ascending node 


(2), argument of periapsis (@) and the time of periapsis passage (T). See Figure 9. 
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Figure 9. Three-dimensional Coordinates from Bate et. al. (1971) 


Ze Two-dimensional 

In the perifocal coordinate system, oa the semimajor axis (a), eccentricity (e), 
and true anomaly (v) are required to fix a position on an orbit. (Note the argument of 
periapsis is defined to be zero in this coordinate frame.). Other parameters of concern are 
the perihelion distance (q) and the aphelion distance (Q). Figure 10 displays the two- 


dimensional perifocal coordinate system. 





Figure 10. Perifocal Coordinates 
D. KEPLER’S SECOND LAW 
Motion of an object about the primary focus in an elliptical orbit is governed by 


Kepler’s Second Law: the line joining the planet to the Sun sweeps out equal areas in 


equal times. This relation increases the difficulty in the problem solution in that there 


Zi 


exists a trancendental relationship between time and position 


Kepler’s Equation and takes the form: 


E-esinE 
t = —__—__, 
n 
e+cosv 
where cos E = ——————_. 
1+ecosv 


I 


. This relationship is called 


(1) 


(2) 


Va PROBLEM FORMULATION AND SOLUTION 


A. STATEMENT 

Given an impending asteroid impact with the Earth, one would like to know if 
there is an optimum point on the asteroid’s orbit where an impulse may be applied to yield 
the greatest achievable separation distance at the closest point of approach for a fixed 
impulse. Additionally, it 1s desirable to determine if is there an optimum direction 
associated with the applied impulse that further increases the separation distance. 
B. ASSUMPTIONS 

In order to proceed with an approximate solution method, simplifying assumptions 
were made. The first assumption was that of two-body motion. That is, the Sun is the 
mass center about which the asteroid and the Earth orbit. Furthermore, the asteroid and 
Earth do not gravitationally interact with each other. It was also assumed that there were 
no external perturbing effects on the orbits due to non-gravitational forces other than the 
applied perturbing impulse. All orbits were assumed to be coplanar, which yields a two- 
dimensional problem. The perturbing impulse 1s assumed to occur instantaneously. The 
asteroid is assumed to be one of the near Earth objects and hence in an elliptical orbit 
around the Sun. Hyperbolic and parabolic orbits were not considered for this analysis. 
Finally, it was assumed that the Earth 1s in a perfectly circular orbit at LAU. 
C. TEMPORAL CONSIDERATION 

In determining the separation of a NEO from the Earth, trme becomes the 
dominant factor in solving the problem. The relative phase of each of the orbiting bodies 
determines whether the bodies will collide. Hence, the Earth-to-NEO separation distance 
is the quantity of concern. Changing the orbital elements becomes secondary to changing 
the asteroid's orbital phase with respect to the Earth. 
D. METHOD OF SOLUTION 

Solution to the above problem may be achieved by use of a numerical simulation 
scheme where the orbital equations of motion are integrated from some initial condition 
forward in time. However, such a simulation is time consuming (on the order of several 
minutes per impact scenario) and therefore limits the scope of any analysis. The 


assumptions stated above allow for use of analytic elliptical orbit equations. Building a 
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solution based on these analytic two-body equations offers a greatly reduced time for 
solution (on the order of a couple seconds per impact scenario) and therefore a greater 
scope of analysis may be performed. A Mathworks MATLAB model was constructed to 
numerically execute the following method. A numerical integration simulation was also 
developed in order to validate the analytic method. The solution description below is 
quite general. For a detailed solution description see Appendix A. The MATLAB script 


that corresponds to the solution method is displayed in Appendix B. 


iL Geometry 

A general description of an elliptical orbit intersecting a circular orbit is used in the 
problem solution. This description suffices for all planar mtercept scenarios. A rotation of 
the perifocal coordinates may be required to bring the model in alignment with the 


particular problem, but the relative geometry remains unchanged. Figure 11 demonstrates 


Q 


Figure 11. Equivalent Impact Scenarios 


this equivalence. 


To uniquely fix an intercept scenario with the above assumptions, only the impact 
true anomaly and orbital eccentricity need to be defined. In the case of planar orbits, the 
perihelion direction, as defined by original asteroid elliptical orbit, defines the principal axis 
from which the impact true anomaly is measured. Implicit in the above impact location 


description is the assumption of an Earth orbital radius of 1 AU. 
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ee Solution Flowchart 


Figure 12 depicts the steps in the problem solution method. 


Specify: 
True Anomaly at Impact, 
Asteroid Eccentricity, 
Time from Impact of Impulse. 


Determme: 
Unperturbed asteroid orbital elements. 









Determine: 
Impact time from perihelion with respect 
to unperturbed asteroid orbit. 


Determine: 
Time from perihelion, true anomaly, and 
V for the point of impulse with respect to 
the unperturbed orbtt. 


Apply impulsive perturbation AV 


Determme: 
Perturbed orbital elements from 
T and V+AV. 















Determine 
Time from perihelion and true anomaly 
for the pomt of impulse with respect 
to the perturbed orbit. 


Create a mesh of orbital positions equally 
spaced in time about the mpact pomt on 
the unperturbed orbit. 


Map the unperturbed mesh onto the 
perturbed orbit. 


Map the unperturbed mesh onto the 
Earth’s orbit. 





Determine: 
Earth to asteroid separation distance for 
the perturbed and unperturbed orbits. 


Evaluate the mimmum separation distance 
for the perturbed orbit. 





Figure 12. Solution Flowchart 
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S: Given Conditions 

For a given impact scenario, it is assumed the impact true anomaly, asteroid orbital 
eccentricity are known and the desired impulse time prior to impact is specified. It should 
be noted here that specifying impact true anomaly and orbital eccentricity fixes the semi- 
major axis and hence, the orbital period. If either parameter is changed, the semi-major 
axis changes. That is, as impact true anomaly increases from 0 to 7 the semi-major axis 


decreases, thus decreasing the orbital period. 


4. Unperturbed Orbital Elements 
From the given conditions and the known Earth orbit the perihelion distance, semi- 


major axis, and orbital period for the unperturbed orbit may be determined by: 


_ (1 + eCOS v,) 


P l+e 


2 


* 


= =e 


2 
and P=. 
n 


where n = ee 


The time from perihelion of the impact with respect to the unperturbed orbit may 





a 


5. Impact Condition 


be determined by Kepler’s Equation as given in Equations (1) and (2). 


6. Initial Impulse Condition 
The time from perihelion of the impulse with respect to the unperturbed orbit may 
be determined by Figure 13 and the equation: 


(Impulse Time)perinetion = (Impact Time)perihetion - (Impulse Time) tmpact Time- 
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Figure 13. Impulse to Impact Time Relation 


The true anomaly at impulse may be determined by inverting Kepler’s Equation 
such that true anomaly becomes a function of time. This approach requires an iteration on 
eccentric anomaly. 

The impulse position, r, and velocity, V , may now be found from: 

a(1 —e is 


————— (3) 


1+ecosv’ 


ft =(rcos v)x +(rsin v)/, (4) 
and v7 = Ele sin v)X + (e + cos v)y| : 


where p = a(1 ~ e”). 


Ds Perturbation 

An orbital perturbation, AV, is applied to the asteroid at the impulse position Tr. 
This yields a perturbed orbital velocity V + Av. 

8. Perturbed Orbital Elements 

From the impulse position r and perturbed velocity V + AV the perturbed orbital 
elements may be determined from: 


a 


h=frx(v¥+Av), 


a 


( V+ sop) —{F-(V+ AvV)} (V+ 10) ' 


and Equation (4). 


a 


a. Perturbed Impulse Condition 

Substituting the perturbed orbital elements into Equation (3) enables the true 
anomaly at impulse to be determined with respect to the perturbed orbit. The time from 
perihelion of impulse with respect to the perturbed orbit may be found from Kepler’s 


Equation. 


10. Orbital Positions at Impact Time 

To determine the separation distance of the perturbed asteroid from the Earth, the 
positions, and times, of the unperturbed asteroid orbit about the impact point must be 
mapped onto the corresponding positions and times of the perturbed asteroid orbit. That 
is, a one-for-one correlation between where the asteroid would have been if not for the 
perturbation and where the asteroid is due to the perturbation must be developed. This is 
the key step in determining the effect of the impulse perturbation. 

The true anomalies and time from perihelion of the perturbed and unperturbed 
orbits are related by Kepler’s Equation and not a simple function. The approach used to 


achieve the required mapping is as follows: 


a) Unperturbed Orbital Positions 

Develop an evenly spaced window, or mesh, of time about the impact 
position wide enough to include the perturbed orbits’ minimum Earth separation. (A first 
estimate of the required width of this mesh is achieved by a numerical simulation. From 
the numerical simulation and model development an interval width of + 1.5 x Av x 
Impulse Timeimpact Time WaS uSed to provide a sufficiently wide mesh to obtain a solution 
without excessive computation time.) By inverting Kepler’s Equation, the true anomalies 
of the mesh points may be determined. From the true anomalies of the mesh points the 


orbital positions may be found from Equations (3) and (4). 


b) Perturbed Orbital Positions 

From the relationship of the time of impact (known) with the time of 
impulse (also known) the center of the mesh may be determined for the perturbed orbit. 
Again, the true anomalies of the mesh points and the orbital positions for the perturbed 
orbit may be determined by use of Kepler’s Equation and Equations (3) and (4). This now 
yields the asteroid orbital positions due to the perturbation. 
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Cc) Earth Orbital Positions 
In the same fashion, the Earth orbital positions corresponding to the mesh 
points may be determined. However, since a circular Earth orbit was chosen, it is easier to 


relate the mesh times to orbital positions by use of the Earth’s orbital mean motion. 


11. Earth to Asteroid Separation Distance 

With the perturbed asteroid’s orbital positions and the corresponding Earth orbital 
positions known, it is a simple matter to determine the Earth to asteroid separation 
distance at each of the mesh points from: 

Ry ep = ba = x.e) a (Y 5° Si 

12. Minimum Separation 

From the above orbital separation distances, the minimum may be determined. It is 
useful to express this separation in terms of Earth radii. This enables a quick evaluation of 
whether a sufficient separation was achieved to cause a “miss”. As with the analysis 
performed by Ahrens and Harris (1994) the resultant separation distances scale linearly 
with the applied impulse magnitude. 


13. Sample Model Output 

A sample of the MATLAB model output is shown in Figures 14 and 15. The 
impact scenario is for the case of an impact occurring at a true anomaly of 30°, an asteroid 
orbital eccentricity of 4, and a perturbing impulse of 1 ms’ in the direction of the velocity 
vector occurring 0.47 asteroid orbits prior to impact. Figure 14 displays the initial 
conditions of the scenario at the time of impulse. The impact position and the asteroid and 
Earth positions at the time of impulse are also displayed. Times prior to impact are 


displayed on the asteroid orbit in tenths of an orbital period. 
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Impact True Anomaly = 30 Orbital Eccentricity = 0.5 


Semi-major as = 1.91 AU 
dw = 1 mis 
dvn = O més 


Y (AU) 


® -0.47 orbits 





Figure 14. Initial Conditions 


Figure 15 shows the effect of the impulse on the asteroid near the impact point. 
The unperturbed asteroid position is shown achieving impact conditions at zero Earth 
radii. The perturbed asteroid trajectory is shown approaching an impact condition, but 
instead reaches a minimum separation (indicated by the 'x') and then increases in 


separation. 


Perturbed v. Unperturbed Asteroid to Earth Separation 
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Figure 15. Perturbed Asteroid to Earth Separation 
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The minimum separation indicated in Figure 15 is the final output of the model. 
The figures were generated from the model working variables. This enables the model to 
be incorporated into a controlling routine allowing for repeated simulations that sweep 


over ranges of the input conditions. 
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VI. PROGRAM VALIDATION 


Prior to making an analysis of various impact scenarios, the solution method and 
MATLAB model had to be validated. This validation was achieved by use of a simple 
numerical integration simulation and by comparison to approximate analytic solutions of 
nearly circular orbits. 

A. NUMERICAL SIMULATION 

The numerical integration simulation was developed usingthe Mathworks 
MATLAB and SIMULINK numerical processors. The SIMULINK diagram and 
MATLAB script files that support the SIMULINK model are displayed in Appendix C. 

The two-body equation of motion integrated by the model is: 


tes 
re 


Expressed as a two-dimensional system of linear differential equations in perifocal 


coordinate form: 


Numerous orbital scenarios were run to ensure the numerical integration was 
performing correctly. A fourth order Runge-Kutta integration scheme was used for the 
simulation acting on a second order differential equation. This yields a solution that is 
analytically exact and accurate to the numerical precision of the computer microprocessor. 
For the present analysis, the computer microprocessor was a 100 MHz, 32 bit, Intel 
Pentium. For the cases chosen for orbital modeling verification, integrating around one 
orbit resulted in an ending position the same as the starting position with a relative error of 


2x10'° (3x10” km error at 1.5x10° km, 1 AU). 
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Once the numerical integration orbit models were validated, numerous perturbation 
simulations were made to provide a reference data base for the analytic method. It was 
found that the analytic method numerical model was in excellent agreement with the 
numerical simulations. Differences arose only from the difference in mesh size between 
the two solution methods, with the analytic solution having a finer mesh interval. 

B. APPROXIMATE ANALYTIC SOLUTIONS 

For circular orbits, Ahrens and Harris (1994) performed an approximate analysis 
that yields expressions for maximum orbital separation due to impulses applied both 
normal to and parallel to the orbital velocity vector. For the case where the impulse is 
applied normal to velocity vector, the maximum separation is found half an orbital period 


later with a magnitude 6... = 2AvP / 2. This case is shown in Figure 16. 


ea Snax * 2ZAVP/n 


Perturbed 
orbit .” 





Vv 
Figure 16. Impulse Normal to Velocity Vector 


For the case where the impulse is applied parallel to velocity vector the separation 
a full orbital period later is found to be 6 = 3AvP per orbit. That is, for this case the 
separation increases by 6 for every orbit after the single impulse. This case is shown in 


Figure 17. 
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Figure 17. Impulse Parallel to Velocity Vector 


Evaluation of these scenarios by numerical simulation verifies the above 
approximate solution. The separations between the perturbed and unperturbed orbits at 
the described locations does indeed agree very well with the approximate analytic solution. 
Further, evaluation of the above scenarios by the analytic method numerical model using 
an asteroid orbital eccentricity of 10° (nearly circular) and an impulse of 1 ms" yields 
virtually the same result as the approximate analytic solution. Table 2 summarizes the 


performance of the three solution methods for the circular orbit case. 






Approximate Analytic Numerical | Numerical 









Impulse Direction | Analytic Solution Method Simulation 
Normal to Vv 3AVP = 14.8 M48 8 14. M48 


(Separation in Earth Radi) 
Table 2. Validation of Solution Method 
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VII. ANALYSIS AND RESULTS 
A. ANALYSIS 


The preceding analytic numerical method for solution of the Earth to asteroid 
separation problem was incorporated as a function into a routine that sweeps over impulse 
direction and time of impulse. In this manner hundreds of solutions may be generated in a 
few minutes. Appendix D contains numerous analyses for impact scenarios that may 


occur around the Earth’s orbit. 


i Long Time Response Behavior 

For a perturbing impulse time well prior to impact the observed behavior resulting 
from the orbital dynamics is not linear. The Earth-to-asteroid separation achieved by an 
impulse is strongly dependent on the location of the impulse on the orbit as well as the 
direction of the impulse with respect to the orbital velocity. The typical behavior of an 
impulse is depicted in Figure 18. Each point on this plot represents the minimum 
separation point depicted on Figure 15. In Figure 18 the vertical axis represents the 
minimum Earth-to-asteroid separation (in Earth radi) that results from an impulse of 1 m 
s'. The two horizontal axes represent the direction of the impulse and the time prior to 
impact of the impulse. The axis representing direction of impulse is measured from 0° to 
360° with respect to the forward direction of the velocity vector. The sense of the 
direction is that impulse directions between 0° and 180° point inward to the orbit and 
directions from 180° to 360° point outward from the orbit. The remaining axis represents 
impulse time prior to impact as a fraction of the asteroid's original orbital period (e.g., t/P 


= (0.5 corresponds to an impulse one half orbit prior to impact). 
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Impact True Anomaly = 30 deg 
Orbital Eccentricity = 0.5 
dv=1 m/s 
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Figure 18. Earth to Asteroid Separation 
2 Relative Maxima 


Impulse times ranging from 0 to 1’% asteroid orbits prior to impact yield two 
maxima points as shown in Figure 18. Inspection of these maxima points reveal that they 
occur at the time of perihelion passage for the asteroid approximately one orbit prior to 
impact for impulse directions parallel and anti-parallel to the orbital velocity vector. This 


relationship holds true for all cases considered. 


5: Non-perihelion Impulse Direction 

Closer inspection of Figure 18 reveals that for impulse times less than that 
corresponding to one orbit prior to impact, the maximum separation is achieved if the 
impulse direction is not aligned either parallel or anti-parallel to the velocity vector. This 


effect is better shown in Figure 19, which is a contour plot of Figure 18. 
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Figure 19. Contour Plot of Figure 18. 


A distinct shift in the direction of impulse for maximum separation as the time of 
impulse becomes closer to impact can be seen. This shift in direction arises from the two 
ways in which to achieve the desired separation. If sufficient time prior time prior to 
impact exists, changing the speed of the asteroid on its orbit will shift the phase of the 
asteroid with respect to the Earth and avoid the impact. Ifthe time prior to impact is 
short, changing the direction of the asteroid laterally with respect to its approach to the 
Earth is necessary to avoid the impact. For times prior to impact between these the 
former and the latter, a tradeoff between changing the orbital speed and displacing the 
asteroid laterally on its orbital path exists and is the cause of the shift in the optimal 


impulse direction. 


4. Periodic Growth 


Thus far, only impulse times within approximately one asteroid orbital period of 


impact have been considered. Extending the analysis to beyond this time period 
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demonstrates the manner in which the separation grows as a function of impulse time. 
Figure 20 shows this behavior for the same impact scenario as discussed above. Of note is 
the periodic growth over time and the peak displacements occurring at each perihelion 
point. This particular figure represents a slice of Figure 18 along the 0° impulse direction 
extended from 0 to 10 orbits prior to impact. 
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Figure 20. Periodic Growth of Separation Distance 
B. RESULTS 
The collection of the numerous impact scenarios modeled by the above method 
may be found in Appendix D. A study of these numerous model results yield the following 


general conclusion. 


1. Optimum Impulse Condition 

Assuming the optimum impulse condition is achievable in terms of a “real”? mission 
sense (that is, the booster technology and energy delivery mechanism exist for asteroid 
mitigation), the optimum impulse point is located at the perihelion of the original asteroid 


orbit at least 2 orbit prior to impact. If time prior to impact permits, impulse at perihelia 
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multiple orbits prior have an even greater effect on separating the asteroid from the Earth 
at the given impact time. However, impact prediction becomes a problem in such cases as 


the validity of orbital prediction models in a general n-body problem comes into question. 


Z. Optimum Impulse Direction 

If, due to the late time of detection, it is not possible to achieve even the first 
relative maximum for the optimum impulse condition, then there exists an optimum 
umpulse direction at the time of impulse that maximizes the Earth to object separation. For 
fractions of an orbit from 0.2 < (t/P) < 0.9 there is a shift in the direction of impulse 
toward the orbital inward normal (for impulses that increase the orbital speed) and toward 
the orbital outward normal (for impulses that decrease the orbital speed) that yields a 


maximum in separation achievable. 


3. Time of Arrival Consideration 

Given the periodic nature of the impact problem, there exist conditions where it 
may be beneficial to delay a deflecting impulse until an optimal impulse condition occurs. 
For the scenario corresponding to Figure 20, if time permits delivery of an impulse two 
and one half asteroid orbits prior to impact, it would prove more advantageous to delay 
the impulse until only two orbits prior to impact in order to maximize the effect of the 


impulse. 


4. Detection Consideration 

The difficulty in realizing the use of an optimal impulse condition is that the 
detection ofa colliding object may occur too late to achieve the most desirable condition. 
The earlier the detection the better the chance of deflection with a much smaller imparted 
energy. Unfortunately, the current search programs and record of detection have been 
yielding very short response times for NEO's having very close approaches. The current 
range of times for detection has been on the order of five days prior to closest approach to 


detection only after Earth passage. 
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VIII. APPLICATION TO A REAL CASE 


A. IMPULSE ACHIEVABLE 

The analysis of the energy coupling between an explosive yield and an asteroid 
performed by Ahrens and Harris (1994) shows that it may be feasible to deflect a NEO 
with an impulse magnitude from 1 cms’ to a few ms’ for a globally threatening object of 
about 1 km in diameter. This impulse may be achieved using one of several methods. 
However, due to the relatively short warning times considered in this thesis, a nuclear blast 
appears to be the most efficient energy delivery mechanism. The blast may be a standoff 
detonation or surface detonation. 

For the standoff detonation, the required explosive yield, W, may be determined by 
the approximate expression: 

10° AvD? 
Sr 

from Ahrens and Harris (1994), which has been modified to express yield in kT of 
equivalent TNT. The impulse, Av, is expressed in ms” and the asteroid diameter, D, is in 
km. The efficiency of neutron production, n, from the nuclear blast lies between 0.03 and 
0.3. The standoff blast efficiency factor, A, is taken for an optimum standoff distance of 
0.4 asteroid radii with an efficiency of approximately 0.3. An order of magnitude analysis 


of the above oot is presented in Table 3 
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Table 3. Impulse and Diameter v. Standoff Explosive Yield 










For the surface detonation, the required explosive yield, W, may be determined 
from the approximate expression: 


W = 4x10” AvM,,.0- 
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also from Ahrens and Harris (1994) which again has been modified to express yield in kT 
of equivalent TNT. The impulse, Av, is expressed in ms” and the asteroid mass, Mygo, is 


in kg. An order of magnitude analysis of the above approximation is presented in Table 4. 


Impulse (ms) | 0.1 km diameter 10 km diameter 
600 T 600 kT 600 MT 


60 kT 60 MT 60 GT 


Table 4. Impulse and Diameter v. Surface Explosive Yield 













B. TOUTATIS 

Asteroid 4179 Toutatis will have multiple close approaches with the Earth over the 
next 15 years. It is of interest to apply the above analysis and methodology to Toutatis as 
if it were going to impact the Earth. 

To perform this analysis it must be assumed that the orbit of Toutatis lies in the 
ecliptic plane. This is not far from the true geometry of Toutatis' orbit where the orbital 
inclination is 0.47° out of the ecliptic plane. From the catalog of NEA orbital elements 
listed compiled by Tholen (1995), the semi-major axis and eccentricity of Toutatis are 
currently listed as 2.5154 and 0.6361, respectively. Assuming that Toutatis and the Earth 
will collide and that the Earth is in a circular orbit at 1 AU yields an impact at +38.53° 
with respect to the perihelion of Toutatis. This is determined by solving Equation (3) for 
true anomaly. Figure 21 depicts this relative orientation of Toutatis with respect to the 
Earth. For the following analysis, the +38.53° impact location is chosen for modeling 
purposes. 
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Impact True Anomaly = 38.53 deg, Orbital Eccentricity = 0.6361 


Semunajor axis = 2.52 AU 


Y (AU) 





5 4 3 o a 0 1 
X (AU) 


Figure 21. Relative Orientation of Toutatis Impact 


Modeling the Toutatis collision in the manner described above yields the results 
displayed in Figure 22. From the JPL and NASA Photo Caption (1993) and Press Release 
(1996), the size of Toutatis is estimated to be approximately that of two attached spheres 
having diameters of about 4 km and 2.5 km. If both masses are combined, the effective 
spherical diameter is approximately 4.3 km. Using a mean asteroid density of 3000 kg m” 
results in a mass for Toutatis near 1.25x10'* kg. From the previous impulse analysis, an 
explosive yield of about 5 MT 1s required in the case of a surface detonation and an 
explosive yield from 9 to 90 MT, depending on neutron production, is needed for a 
standoff detonation to achieve a 1 cms" change in orbital speed. 

Using the 1 cms" orbital speed change determined above, a maximum separation 
of 1.64 Earth radii is the result of an impulse delivered 1.02 orbits prior to impact 
(perihelion passage of the prior orbit). While the separation 1s sufficient to cause the 
Toutatis to miss the Earth in this scenario, a greater margin of safety would be desirable. 


A larger separation may be achieved by either delivering an impulse greater than 1 cm s"' 
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or delivering the impulse at an earlier perihelion passage. Recognizing the extremely large 
explosive yield requirements for increased impulse magnitudes, it would appear preferable 
to deliver the impulse at an earlier time. This type of analysis demonstrates the necessity 


for detection of threatening asteroids many orbits prior to impact. 


Impact True Anomaly = 38.53 deg a 
Orbital Eccentricity = 0.6361 a : 
dv = 0.01 mis ee 


Separation (Earth Radii) 


360 





Figure 22. Impulse Effects on Toutatis-Earth Impact 
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IX. FURTHER CONSIDERATIONS 


The above analysis shows promise as a tool for rapidly evaluating numerous 
scenarios for deflecting an asteroid that is going to impact the Earth. The possibility for 
further investigation utilizing this method presents itself in analyzing longer response times 
prior to impact, generalizing the method to a three-dimensional method and eccentricities 
other than “2. Additionally, more work is needed in analysis of very short response time 
impulse effects. 

A. LONG RESPONSE TIME 

The analysis presented above and in Appendix D have been performed primarily 
for impulse times between 0 and 1% orbits prior to impact. A few models were made for 
one impulse direction at times ranging from 0 to 10 orbits prior to impact. However, this 
investigation needs to be pursued further in search of general trends other than maximum 
separations occurring at perihelion points. 

B. THREE-DIMENSIONAL ANALYSIS 

The current model and method apply only to two-dimensional scenarios. It is of 
interest and merit to further generalize the analysis to the three-dimensional case. This will 
allow for orbital inclinations out of the ecliptic plane and better simulate a variety of real 
scenarios. 

OF ECCENTRICITY VARIATIONS 

In the above analysis, other than for Toutatis, the orbital eccentricities have been 
maintained at %. An investigation into the effects of more circular orbits and highly 
eccentric orbits is in order. 

D. SHORT RESPONSE TIME 

The method presented is derived from a two-body representation of a more 
complicated physical system. This method is not valid for very short response times when 
the impacting object is within the Earth's sphere of influence. A further investigation is 
desirable to assess the effects of impulses in the three-body problem that arises when the 


object is detected close to the Earth. 
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APPENDIX A. DETAILED SOLUTION METHOD 


Given an Earth-asteroid collision with the following properties: 


A circular Earth orbit with semi-major axis (radius), a, = LAU 
(orbital eccentricity, e, = 0); 
an asteroid orbit with orbital eccentricity €,, ; 


the true anomaly at time of impact, v impact,U = V inpact,E > 


t 
and the time of the perturbing impulse, AV, prior to impact & 


U Av 


Find the minimum separation of the asteroid from the Earth in the vicinity of the old impact point. 


55, 











Variables: 


Nomenclature 


semi-major axis 

eccentric anomaly 

orbital eccentricity 

specific angular momentum 

fraction 

number 

orbital mean motion 

orbital period 

orbital parameter 

asteroid to Earth separation distance, km 
orbit to Sun radial distance 

velocity 

X-position 

y-position 

impulse quantity 

angle w.r.t. x-coordinate axis 
gravitational mass parameter 

true anomaly 

asteroid to Earth separation distance, Earth radii 


shift in periapsis direction 


fraction of orbit from periapsis 


a7. 


Subscripts: 


days 


E 


pos 
Range 
sep 


Sun 


number of days 

Earth 

i” element 

increment 

impact position 
minimum 

orbit count 

perturbed asteroid orbit 
periapsis 

position 

set of all i’s 

separation 

Sun parameter 
unperturbed asteroid orbit 
x-component 
y-component 

impulse 

parallel component 


normal component 


incremental amount 


incremental amount 
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Elliptical Orbit Relations 


Mean motion: 


Ul 
n= ay 
Period: 
27 
a 
rh 


Radial distance to gravitational focus (general): 
a(1 _ e?} 
~ 1l+ecosv 


Periapsis radius 


r,(1+ecos v,} 
r= 
: l+e 


Semi-major axis: 





Eccentric anomaly as a function of true anomaly: 


e+cosv 
cos E = —————_- 
1+ecosv 


True anomaly as a function of eccentric anomaly: 


cosE-—e 
cos v = ———_—__— 
1—ecosE 


Time since periapsis: 


E-esinE 
e n 


True anomaly: 
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COS V = = 


a(1 — e”} 1 
re © 


Position at a point in perifocal coordinates: 
ft =(rcos v)X + (rsin v)y 


Velocity at a point in perifocal coordinates: 


v= JE[c sin v)x + (e + cos vv] 


Parameter of an orbit: 


Specific angular momentum: 


ee, 


h=rxv 


Eccentricity vector in perifocal coordinates: 


s=4|(1r-E}r--9r 
ht Ir] 


Constants 


LAU = 1.4959787x10* km, astronomical unit 

G = 6.67259x10~” km*kg“'s~, gravitational constant 

M,,. = 1.9891x10™ kg , mass of the Sun 

Uoun = GMg,,, = 1.327124399355x10" km’s~, Sun gravitational parameter 


R, = 6.37814x10° km , Earth radius 
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Detailed Solution Method 


I. Determine unperturbed asteroid orbital elements. 


a. Perihelion radius: 


a,(1 + €,, COS ei 


Tr = 
p,U > 
lewe.. 


where Ts mpact,U =a, 


b. Semi-major axis: 





oo, = ae km 
wae, > 
c. Mean motion: 
Ue 
n, = - ,rad s7? 
ay 
d. Period: 
Qn 
Ny 


IJ. Determine the conditions at impact with respect to the unperturbed asteroid orbit. 
a. Eccentric anomaly at impact: 


Cy + COS VinractU 


E. = cos! 
impact,U 
1+e,, cos v es 


(cos” principle values are 0< 89 <7) 


b. Fraction of orbit since (or prior to) perihelion of impact: 


t E snpact,U - ey sin E se mpact,U 
lee. aT 
umpact,p 
t 
if Vee Os Pp. <0 
U impact,p 
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III. Determine the conditions at impulse with respect to the unperturbed orbit. 
a. Fraction of orbit prior to (or since) perihelion of impulse: 


Fg“ “Wed 


Av, impact,p Av ,impact 


b. Eccentric anomaly at impulse (solve via Newton-Raphson iterative method): 


( t Eau 7 Gy COSE,, y 
= = 


P, 27 


Av,p 
c. True anomaly at impulse measured with respect to unperturbed perifocal coordinates: 
cosE 


Vie = COSe [eB 20. au | 
Av U 
i—e,, cose. y 


(cos” principle values are 0 < 0 < 71) 


t 
es a IN rit at Ee 


U Av,p 


i 
if-l< a < “5? then Vayu = 2n(N ore = 1) + Vayu 


r ex 0, 2n 
<F 


if — obit «0, then Va, y =20N 


NO | = 


orbit __ V av,U 


7 Cub 0, 2x 


if 0 S$ te < 2 > then V av,U 2TN awit a Vav,U 
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™ SA On2r 


aI 
if eS Foy <1 theme y= 2n(N ori as 1) — Vav,u 


2 aa or 
7 cD) 0,27 
pit 


d. Distance to focus (Sun) from impulse location: 


ay(1-e%) 


TayU i ? km 
(1 + 1, COS Vav.u) 


e. Velocity components at impulse with respect to unperturbed perifocal coordinates: 


Hoon : et 
een = eens) SUNY av.» Kms 
a,li-ez 
Hoon -] 
Vy avu = fe, + cos V av)» km S 
ay\l-e; 


f. Velocity direction with respect to unperturbed perifocal coordinates 


pre -1 ( Vv ee 
Ws x,Av,U 


(achieved utilizing “atan2(y,x)” numerical routine) 
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IV. Now perturb the asteroid with respect to velocity vector direction. 


a. Impulse velocity components with respect to unperturbed perifocal coordinates: 


Av, = Av, cosa — Av, sina , km s” 
Av, = Av, sina + Av, cosa, km s” 
where, AV, is in the direction of the velocity vector 


and, AV, is normal to the velocity vector in a right hand sense 
b. Velocity components after perturbation with respect to unperturbed perifocal coordinates: 


took -1 
May = YAU 1 AV, ? km s 


a 7 
V,avP = Vyavu + AV: km s 


Av, sine 


Av , COS —» 
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V. Determine perturbed asteroid orbital elements from Fr and V. 


a. Position vector with respect to unperturbed perifocal coordinates: 


TavP 1 tae Ty av,u* ats Ty Av,UY +0z , km 


/ 


where, Ty avu = Tay COS Vayu 
and, Ty Av,U = Pay,U sin Vav,U 


b. Velocity vector with respect to unperturbed perifocal coordinates: 


aA 


— aS A A =| 
V av,P - V x Av,PX 1 Vy av PY +0z , km s 


c. Specific angular momentum vector with respect to unperturbed perifocal coordinates: 


= _ Dat 
h, = Tayp X Vay p» km’s 
d. Perturbed orbit parameter: 
a2 
|h,| 


Pp= We a a,(1-e2), km 





where lho? =h,-h, 


e. Perturbed eccentricity vector with respect to unperturbed perifocal coordinates: 








Ol 
II 


1 U 
P | A ’ | f ( , : e ) s 
V,F [r ; | Av,F Av,P Av,F Av,P 


= 


= 2 = — 
where Boer = Ving PY Ee 


and | | = Taw u 


f. Perturbed orbital eccentricity: 


g. Perturbed orbit semi-major axis: 


Pp 


=e km 
ee 


ap 


h. Perturbed orbit mean motion: 
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Us 
Dc | ae TAdise 
ap 


i. Perturbed orbital period: 


VI. Determine the angle between the perturbed and unperturbed orbital eccentricity vectors. 
a. Unperturbed orbital eccentricity vector: 
éy =eyx+ Oy +0z 


b. Rotation of eccentricity vector due to impulse: 


e.-é 
av = cos( : | 
Cply 
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VII. Determine the impulse conditions with respect to the perturbed orbit. 


a. True anomaly at impulse with respect to the perturbed orbit: 


a,(1-e3) 1, | 


-1 
Vv = cos7?| ——————- — 
2) r. .e e 


Av,P ~P P 


(cos” principle values are 0 < 8 < 77) 


1 
ifr ko ee < 9? then V av,P a 2n(N ocvit a 1)+ Vav,P 


<F 


orbit 


i 


Nl eR 


< 0 9 then V av,P = 2nN orbit V av P 


0, 2x 


\ 


1 
if O< Facey < 2 ’ then V AvP = 2nN orbit ae V av,P 


E 0, 2x 


() 
y 


1 
if 5 < Fo <1, then Vy, p = 2n(N Deets 1) Vinee 


O72 


é 
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b. Eccentric anomaly at impulse with respect to perturbed orbit: 


ee €, +COSv, 5 
Be 1+e, cos v 
P Av,P 
(cos principle values are 0 < 8 < 71) 


c. Fraction of orbit of prior to (or since) perihelion of impulse with respect to perturbed orbit: 


(+ _ Eee ctu ee 


Pe Rs Qn 
1 t t 
if -1<F_,, <——x, then (. = (Noize —1)+ ft 
2 P, Av,p P~ av p 
1 t t 
if 5 F pit < O, the =) Ne ( 


Pye t 
if O< <F, ees >? , then ae 


Pp 


Av,p 





tt t t 
if 5 SF <1, then Pp) = (N jae + 1)- a 


Av,p P@ Av,p 
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VIII. Determine a range of times, and hence positions, of the unperturbed asteroid in the vicinity of the 
impact point on the unperturbed orbit for mapping to corresponding times, and positions, of the perturbed 
asteroid on the perturbed orbit. 

a. Set time period: 


N = =) sdays 


days 
b. Set number of positions to map during time period: 
n pos = 201 
c. Establish limits of time interval: 
sf t [Son 2Abes) a 
Pe - lee lday Thr 
d. Number of increments in time period: 
Nn imc es Nl nos _ i 


e. Incremental time step: 


te) les 


f. Now have Ngo positions from -Ngzys to +Naays in steps of 2Naay/Nine parts of a day represented 





as fractions of a complete orbit. 


Ge), Am) ag) sa) mama 
ay =—Al— | = ||] S ——], instepso = 
Pe ee ee P,, P,, 


g. Center this range of positions on the fraction of orbit since (or prior to) perihelion of impact: 
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% 


h. Eccentric anomaly at each time coordinate (solve via Newton-Raphson iterative method): 


( t Eiu au cosE; y 
- 27 


(a), 


Range,impact Range impact,p 


le 


U i,Range,impact 


i. Determine the true anomalies at each time coordinate with respect to the unperturbed perifocal 


coordinates. 


ee Conte aiea 
Lo = 
=e), Conese 


(cos principle values are 0 < 9 < 72) 


t 1 
; Fe ee 
1,Range, impact 
t 
if Pp < 0 \thene yey 
U i,Range,impact 
art 1 
if Pp. > >? then V;y =2N—-—Viy 


U i,Range impact 
j. Determine the radial distance to the focus (Sun) for these time coordinates: 


a,(1—e?,] 
ae 
eer, COSY, 


k. Determine the distance components at each time coordinate with respect to the 


unperturbed perifocal coordinates: 


Xiu =T,y COS V; y, km 


1, 


Yiu =ySMvyy. km 


1 
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IX. We now need to determine the true anomaly and fraction of orbit on the perturbed orbit that 
corresponds to the true anomaly and fraction of orbit at the impact position on the unperturbed orbit. 
Determine the positions and orbit fractions on the perturbed orbit that correspond to the Npos 
positions and orbit fractions of the unperturbed orbit. 


a. Fraction of orbit since (or prior to) perihelion of the “impact” position with respect to the 


le) 


Av 


perturbed orbit: 


a), 


ct Av,p 
b. Establish the n,,; orbital positions and fractions of orbit on the perturbed orbit that correspond 
exactly to the n,,, positions of the unperturbed orbit: 
f t RP t 
oS) ee 
P Range,impact U Range P P impact,p 


c. Eccentric anomaly at each time coordinate (solve via Newton-Raphson iterative method): 


t E,p—e, cosE,;, 
a 27 
i,Range,impact 
d. Determine the true anomalies at each time coordinate with respect to the perturbed perifocal 
coordinates. 


cosE,,—€p | 


-j 
V., =CcOSs ( 
ie 


(cos principle values are OS 9 S71) 


t 1 
if | — <—-—,then v,p =—2N+V;p 
Pe ? , , 
i,Range,impact 
Lainie 
if Pp. <0, then V;p =—V;p 


P i,Range,impact 
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if es > =: then V;p = 2% —Vip 
P~ i Range,impact 
e. Shift true anomalies from perturbed coordinates to unperturbed coordinates: 
Vipu = Vip + Av 
e. Determine the radial distance to the focus (Sun) for these time coordinates: 


_ ay(d aoe) se 


Ww 14+e,c0Sv,_p” 


f. Determine the distance components at each time coordinate with respect to the 
unperturbed perifocal coordinates: 


Xip =Tip COS ViPu> km 


Yip =T,p SMV; py. km 


X. Now we need the positions of the Earth at the same n,,, positions surrounding the impact position. 


a. True anomaly of the Earth at each of the n,,, positions: 


t 
VRange,E = V impact,E fs & (rene) 
U Range 
b. Determine the distance components of the Earth at each time coordinate with respect 


to the unperturbed perifocal coordinates: 


Xip =a, COSV,,,km 


Viz = sin V; > km 
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t t 
XI. Now must handle the special case of | < of | ; 
UZ ay 18 


For this case, set: 


Xip =Xiv 

Dear ) 10 

while | —— Ae 
Sar Av Pa 


XII. Define the separation distances of the perturbed and unperturbed asteroid orbital positions from the 
Earth orbital positions. 
a. Perturbed asteroid to Earth separation distance: 


2 2 
Roop = oe —x; «| +(¥ip “Yiz) » km 


b. Unperturbed asteroid to Earth separation distance: 


2 2 
Rou = Ri. -x;.) +(Y.u -Yiz) » km 


c. Scale the separation distances by the.radius of the Earth: 


sep,P 








Prop P = , Earth radii 
E Ke 
wu 
eine = , Earth radii 
E 
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XIII. Finally, select the minimum value over the orbital position range as the minimum separation 
distance. 


a. Minimum perturbed asteroid to Earth separation: 
(scp? = min(p...p). Earth radii 


b. Minimum unperturbed asteroid to Earth separation (By the definition of the problem this is 


zero, but this provides a good check of the method.): 


(p.. ss) = min(p..u , Earth radi 
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APPENDIX B. ANALYTIC SOLUTION METHOD, MATLAB MODEL 


WVWVVVVVWVWVV%VVV0%0/0%0 0% 00% VV %%0%%%%MU%%%%%%%%%%%M%%%%M%% 
% Arbitrary Elliptical Orbit Intersecting Circular Earth Orbit 

% 

% MATLAB Script File Name: nwarborb.m 

% 

% Author: LT Jeffrey T. Elder 

% Naval Postgraduate School 

% November 1996 

% 

RWVVV%V%%0%0%0 0% 0% VV MVVO%V%%%%%%%V%V%VVMUVVM%%MUVUV%M%U%%%% 


RWWVWVWVVWVVV0%% VV %0%0 VV Vo VV %%%V0V%% 0% %%%0%% VOV%%%%MUM%%%M%%%VO% 
% Initialize computational workspace 
clear 


YM %%%%%%%%%%%%%%%%%%%%%%%%M%MMMM%M%%MM%MM%MM%M%M%M%M%%M% 
% Define constants 


AU = 1.4959787e8; % astronomical unit 
musun = 1.327124399355el11; % gravitational parameter for the Sun 
epsilon = Se-6; % Newton-Raphson tolerance 


NMUMVVWVVW%%V% VV %V%% 00% VoVo VOVOV V0O%VO%% VV V%VV%VV%V%%%%%%0% 
% Gather input data from user 

nuimpactU = input(‘Enter True Anomaly for Impact (deg): *); 

eU = input(‘Enter Asteroid eccentricity: '); 

tPUdelv = input(‘Enter (t/P) of Impulse Prior to Impact: "); 

del V = input('Enter the Magnitude of the Impulse dV, (m/s): '); 

thet = input(‘Enter the Direction of the Impulse wrt V, ccw+ (deg): '); 

dvvi = del V*cos(thet*p1/1 80); 

dvni = del V*sin(thet*pi/1 80); 


ARUHUVUV%%UUWM%WUWUVV%W%MUMVUVVVWVW%%%%%%V%%0%%MV%%0%%0 V/V V0 VoVVo%V0Vo/0/0/0 
% Allow for mode! performance evaluation 

Yoflops(0) 

Ytic 


%H%H%%HH%HHM%%%%%%%%M%%%%%%%M%%%%%%%M%%%%%%%M%MV%%%%%%Y%0% 
% Set Earth orbital elements 

aE = 1.0*AU; 

nE = sqrt(musun/aE‘3); 


%™%H%HUH%H%%%%%%M%%%%%%%MMMM%W%M%% VV %%V%V0V0V%V%%VOVoVoVOVOV%V0V0Y%0 
% Convert impact trie anomaly to radians 
nuimpactU = nuimpactU*pi/180; 


%H%%%%HHU%%%%%%%%%%%VMMMUMM0%%0%%0%%%V%%%%0%VoV0 V0 Vo YoYo YoVoV0V0V0/0 V0 
% Determine unperturbed asteroid orbital elements 


rpU = aE*(1+eU*cos(nuimpactU))/(1+eU); % Perihelion radius 
aU = rpU/(1-eU); % Semi-major axis 
nU = sqrt(musun/aU%3); % Mean motion 
PU = 2*py/nU; % Period 
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%Vw%%V%%UW%%%%MW%V%MUUAMMA%MV%VVAMVVV%WVV%V0%0%0%0/oVo/0/o/oV0/0/oVo/o/o/o%o 
% Determine the conditions of impact wrt unperturbed orbit 


EmmpactU =... % ... iS line continuation 
acos((eU+cos(nuimpactU))... 

/(1+eU*cos(nuimpactU))); % Impact eccentric anomaly 
tPUpimpact =... 

(EimpactU-eU* sin(EimpactU))/(2*pi); % t/P of impact wrt perihelion 


if nuimpactU < 0 
tPUpimpact = -tPUpimpact; 
end 


WVMWVAVMMVAHVAVMVAWUWVVMVV%VVAVAVUVAVAVUNVVAVAVVVVVV%% VV VoV%o% 


tPUpimpulse = tPUpimpact-tPUdelv; % t/P of impulse wrt perihelion 
Norbit = fix(tPUpimpulse); % Number of whole orbits 
Forbit = tPUpimpulse - Norbit; % Fraction of orbits 


NVUHV/UwMUHVMUMMMVMAMMAHUAMMMHMMMAMUMUMUVUMMUMMMNMMHMMUMUUVAUAUM%%%% 
% Allow for numerical simulation validation 

tao = tPUpimpulse*PU; % Start time for simulation 

tend = (tPUpimpact+0.125)*PU; % End time for simulation 


HVV%HHMM%%%%MHHMMMM%%%%%MM%MMMMMMMMNMMMMHMNMMMMMMUMMM%MN% 
% Determine true anomaly of impulse 


EimpulseU = pi/4; % Guess eccentric anomaly 

tptemp = Forbit; 

if tPUpimpulse < 0 % allow for +/- tPUpimpulse 
tptemp = -Forbit; 

end 

true = 0; 

while true == 0 % Newton-Raphson iteration 


fprime = 1 - eU*cos(EimpulseU); 

f = EimpulseU - eU*sin(EimpulseV) - 2*pi*tptemp; 
Elast = ErmpulseU; 

EimpulseU = EimpulseU - f/fprime; 

if abs(EimpulseU-Elast) < epsilon, true=1;, end 


end % End N-R iteration 
if tPUpimpulse < 0, EimpulseU = -EimpulseU;, end 
nuimpulseU =... 
acos((cos(EimpulseU)-eU),/... 
(1-eU*cos(EimpulseV))); % True anomaly at impulse 
% If t/P of impulse exceeds 


% one orbit, determine true 
% anomaly for multiple orbits 
if (-1 < Forbit) & (Forbit < -0.5) % Conditions to properly 
nuimpulseU =... % locate true anomaly 
2*pi*(Norbit-1)+nuimpulseU; 
elseif (-0.5 <= Forbit) & (Forbit < 0) 
nuimpulseU = 2*pi*Norbit - nuimpulseU; 
elseif (0 <= Forbit) & (Forbit < 0.5) 
nuimpulseU = 2*pi*Norbit + nuimpulseU; 
elseif (0.5 <= Forbit) & (Forbit < 1) 
nuimpulseU = 2*pi*(Norbit+1) - nuimpulseU; 
end 


16 


LWW WYVWHV/ALDHAD/PA/V/V%%%%7%%VVU%%UW%%%%MV%U%U%HUMHUMUHU%% 
% Determine orbital distances at time of impulse 

rimpulseU = aU*(1-eU%2)/(1+eU*cos(nuimpulseV)); 

ximpulseU = rimpulseU*cos(nuimpulseU); 

yimpulseU = rimpulseU*sin(nuimpulseU); 


HHH %%/HHH%H%M%H%HH%%%%%%%%%%%%H%%H%%H%H%H%%%H%G% 
% Perturb the asteroid 


pU = aU*(1-eU*eU); % Parameter of orbit 
vximpulseU =... 

-sqrt(musun/pU)*sin(nuimpulseV); “ Unperturbed velocities 
vyimpulseU =... 

sqrt(musun/pU)*(eU+cos(nuimpulseV)); 
dvv = dvvi/1000.0; “% Impulse parallel to V 
dvn = dvni/1000.0; % Impulse normal to V 
alpha = atan2(vyimpulseU,vximpulseU); “% Angle of V wrt x-axis 
vximpulseP =... 

vximpulseU+dvv*cos(alpha)-dvn*sin(alpha); % Perturbed velocities 
vyimpulseP =... 


vyimpulseU + dvv*sin(alpha) + dvn*cos(alpha); 


MPAHWDPHHVAVAH/AVV V/V/V VHHVHWMV%AHAMUWHKVHVUHMUMUVUHUMUW*%% 
% Determine perturbed asteroid orbit elements from R and V 


rimpulsePvec = [ximpulseU yimpulseU 0]; % Position vector 
vim pulsePvec = [vximpulseP vyimpulseP 0]; % Velocity vector 
hpvec = cross(rimpulsePvec, vimpulsePvec); “ Angular momentum vec. 
pP = hpvec*hpvec'/musun; % Parameter of orbit 
ePvec = ( (vimpulsePvec*vimpulsePvec'-musun/rimpulseU)*rimpulsePvec... 

- (rimpulsePvec* vimpulsePvec')*vimpulsePvec )/musun; % Eccentricity vector 
eP2 = ePvec*ePvec’; 

eP = sqrt(eP2); % Eccentricity 

aP = pP/(1-eP2); % Semi-major axis 
nP = sqrt(musun/aP‘3); % Mean motion 

PP = 2*pi/nP; % Period 


%%%%%%%H%%U%%H%MMMMMMAMMMMMMMMMMMWMHUHUHMUWVHVMUVUHMWUV%N% 
% Determine angle perturbed orbit perihelion makes wrt unperturbed orbit perihelion 

dnu = acos( (ePvec*[eU 0 0}')/(eP*eV) ); 

if ePvec(2) < 0, dnu = -dnu;, end 


%%%%%%M%%%M%%%%%%H%HHM%M%%MM%%MMM%U MVM AUVHVUWMUWMVAUMVOVO%Y% 
% Determine impulse true anomaly wrt perturbed orbit 
cnuimpulseP = aP*(1-eP2)/(rimpulseU*eP) - 1/eP; 
nuimpulseP = acos(cnuimpulseP); % True anomaly at impulse 
% If t/P of impulse exceeds 
% one orbit, determine true 
% anomaly for multiple orbits 


if (-1 < Forbit) & (Forbit < -0.5) % Conditions to properly 
nuimpulseP =... 
2*pi*(Norbit-1)+nuimpulseP; % locate true anomaly 


elseif (-0.5 <= Forbit) & (Forbit < 0) 
nuimpulseP = 2*pi*Norbit - nuimpulseP; 
elseif (0 <= Forbit) & (Forbit < 0.5) 


1a 


nuimpulseP = 2*pi*Norbit + nuimpulseP; 
elseif (0.5 <= Forbit) & (Forbit < 1) 

nuimpulseP = 2*pi*(Norbit+1) - nuimpulseP; 
end 


%%H%M%H%%%%%%MM%M%%%%MMM%M%M%%%%%%%MVVVV%%0%0%0 YoYo YoYo YoYo %o% 
% Determine impulse time wrt perturbed orbit 
EimpulseP = acos( (eP+cnuimpulseP)/(1+eP*cnuimpulseP) ); 


tPPpimpulse =... 
(EimpulseP-eP*sin(EimpulseP))/2/pi; % t/P of impulse 
% If t/P of impulse exceeds 
% one orbit, determine t/P of 
% impulse for multiple orbits 
if (-1 < Forbit) & (Forbit < -0.5) % Conditions to properly 
tPPpimpulse = (Norbit-1)+tPPpimpulse; % determine t/P of impulse 


elseif (-0.5 <= Forbit) & (Forbit < 0) 
tPPpimpulse = Norbit - tPPpimpulse; 
elseif (0 <= Forbit) & (Forbit < 0.5) 
tPPpimpulse = Norbit + tPPpimpulse; 
elseif (0.5 <= Forbit) & (Forbit < 1) 
tPPpimpulse = (Norbit+1) - tPPpimpulse; 
end 


WWVVWVVVWVV%%UV%V%V%% 70% 0% %0%0%%0V0%0 V0 %0%0% V0 %0 V0 0% 0% 00% %0%0%OVO% 


ndaycoef = 1.5; 

Ndays = ndaycoef*tPUdelv*delV; % npos & Ndays define initial 
% mesh size 

loop = 2; “% Make two passes to ensure 

while(loop >= 1) % ‘mesh’ is properly sized 

npos = max([201 fix(100*Ndays+1)]); % to give accurate minimum 

nuU=[]; % Initialize true anomaly 

nuP=[]; % matrices 


MWAWVVHWWV/0//VVVVVVV% V/V VVVV%%%V% VV %%0%%%0%%%%0%0% Yo 
% Create a mesh of orbit positions about the impact point on the 
% unperturbed orbit 


deltPU = Ndays*24*3600/PU; % Unperturbed mesh half width 
tPURangeimpact =... % Unperturbed mesh properly 
linspace(-deltPU,deltPU,npos)+tPUpimpact;% located in t/P space 
for 1= l:npos % Determine true anomaly of 
true = 0; % unperturbed mesh points 


EU = tPURangeimpact(i)*p1; 
tptemp = tPURangeimpact(i); 
if tPURangeimpact(i) < 0, tptemp = -tPURangeimpact(i);, end 
while true == 0 % Newton-Raphson iteration 
fprime = 1 - eU*cos(EU); 
f= EU - eU*sin(EV) - 2*pi*tptemp; 
Elast = EU; 
EU = EU - f/fprime; 
if abs(EU-Elast) < epsilon, true = 1;, end 
end % End N-R iteration 
if tPURangeimpact(1) < 0, EU = -EU;, end 
nuU(i) =... 
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acos((cos(EU)-eU)./(1-eU*cos(EU))); = % True anomalies of 


if tPURangemmpact(i) < -0.5 
nuU(i) = -2*pi + nuU(i); 
elseif tPURangeimpact(i) < 0 

nuU(i) = -nuU(i); 
elseif tPURangeimpact(i) > 0.5 
nuU(i) = 2*pi - nuU(i); 
end 
end 
rU = aU*(1-eU%2) ./ (1+eU* cos(nuV)); 
xU = rU.*cos(nuU); 
yU = rU.*sin(nuV); 


‘ unperturbed mesh points 
“% Ensure true anomalies are 
% properly located 


% Positions of unperturbed 
% mesh points wrt unperturbed 
% coordinate frame 


PPV VV VVVVY%AHVVWVVHVVLVHVVVVVVUHVMUVUUMUUMVUUVUVUVHUVHUYV%VUH%UUM% 
‘% Map the mesh of orbit positions about the impact point onto the perturbed orbit 


tPPpimpact = tPPpimpulse + tPUdelv*PU/PP; 
deltPP = Ndays*24*3600/PP; 
tPPRangeimpact =... 
linspace(-deltPP,deltPP,npos)+tPPpimpact;% located in t/P space 
for i = 1:npos 
true = 0; 
EP = tPPRangeimpact(i)*pi; 
tptemp = tPPRangeimpact(i); 
if tPPRangeimpact(i) < 0, tptemp = -tPPRangeimpact(i);, end 
while true == 
fprime = | - eP*cos(EP); 
f= EP - eP*sin(EP) - 2*pi*tptemp; 
Elast = EP; 
EP = EP - f/fprime; 
if abs(EP-Elast) < epsilon, true = 1;, en 
end : 
if tPPRangeimpact(i) < 0, EP = -EP;, end 
nuP(i) =... 
acos((cos(EP)-eP)./(1-eP*cos(EP))); 


if tPPRangeimpact(i) < -0.5 
nuP(i) = -2*pi + nuP(i); 
elseif tPPRangeimpact(i) < 0 
nuP(i) = -nuP(i); ? 
elseif tPPRangeimpact(i) > 0.5 
nuP(i) = 2*pi - nuP(i); 
end 
end 
nuPU = nuP + dnu; 


rP = aP*(1-eP2) ./ (1+eP*cos(nuP)); 
xP = rP.*cos(nuPU); 
yP = rP.*sin(nuPU); 


% UP of impact wrt perihelion 
% Perturbed mesh half width 
% Perturbed mesh properly 


% Determine true anomaly 
% of perturbed mesh points 


% Newton-Raphson iteration 


% End N-R iteration 


% True anomalies of perturbed 
% mesh points 

% Ensure true anomalies are 
% properly located 


% Perturbed mesh true 

% anomalies wrt unperturbed 
% coordinate frame 

% Positions of perturbed mesh 
% points wrt unperturbed 

% coordinate frame 


MWH%H%H%%%MH%%M%%%%MM%%%0%%%%M%M%MMM%%1%%MM%M%%%%% 
% Map the mesh of orbit positions about the impact point on the Earth's orbit 


delnuE = Ndays*24*3600*nE; 


WD 


% Earth mesh half width 


nuRangeE =... 


linspace(-delnuE,delnuE,npos)+nuimpactU; % Earth mesh properly located 
% in true anomaly space 

xE = AU*cos(nuRangeE); % Positions of Earth mesh 

yE = AU*sin(nuRangeE); % points wrt unperturbed 


% coordinate frame 


YWHAW*UUYWUWMUVUHUVUUwUWwMMUNMUWWMMUMU%%%UVUMVUVVMUUVVUUVUVVMUMMUMYUAUMUVMVVO%% 
% handle special case of tPUdelv < deltPU 
if tPUdelv < deltPU 
[v,loc] = min(abs(tPPRangeimpact-tPUpimpact+tPUdelv)); 
xP(1:loc) = xU(1:loc); 
yP(1:loc) = yUC1 :loc); 
end 


YHVM%%%%%VMUMMUV%MVUVUWVVM%M%%U%% VV %VVVVVVVV VV VA%V%V%Y% 
% Determine the Earth to asteroid separation in Earth radii 


rhosepU =... 
sqrt((xU-xE).%2 + (yU-yE).%2)/6378.14; % Unperturbed orbit 
rhosepP =... 
sqrt((xP-xE).%2 + (yP-yE).%2)/6378.14; % Perturbed orbit 
[rhosepUmin,inu] = min(rhosepU); % Minimum unperturbed 

% separation (must be zero) 
{rhosepPmin,inp] = min(rhosepP); % Minimum perturbed separation 


WW/WVVWVVVVAV%VV%V%%%%%UVUVW%MUWAMUW%MUUMUMUVAVWVVUUAMUVUMVU%VVUVVV%% 
if Ndays == ndaycoef*tPUdelv*delV % Refine mesh 
Ndays = 1.25*abs(tPURangeimpact(inp)*PU/3600/24 ... 
-tPURangeimpact(inu)*PU/3600/24); 


end 
loop = loop - 1; % Allow for loop exit 
end % End while(test) loop 


MWWWWVHVWVMVMNVHUMWUVMUVMNMNWMMUHMMU%MNMUMM%UMNMMNUMUMANMMMMAV%UMV%%%%V% 
rhosepPmin % Display minimum separation 


PWRwWWVVAWP/VVVWV%V%V%V%VVVVVVMVVVVM%VVMVVUMUWUVV%%VU%VUVUHVWV%VM%%%% 


“toc % Display model performance 
“flops 
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APPENDIX C. NUMERICAL VALIDATION MODEL 


SIMULINK MODEL 
sf ire] 
Clock time 
V1 
TIAB meh 
ry Vi 1 = 
e eapton fs Demux y Unperturbed Orbit 
f 
orb_eomt initial Conditions 
Vx x1 
x yt 
Y Y1 
Demux1 
Pw 
Vio 
TLAB i 
: v2 Perturbed Orbit 
Function a -- 
orb_eom2 Initial Conditions eee 
Vx + dvx X2 
+ 
aoe 
Y Y2 
Demux2 


In order to run the SIMULINK numerical integration model, the analytic method 


model must first be run. The initial conditions for the integrators are specified by the 


working variables in the analytic model. 


81 


MATLAB SCRIPT SUPPORTING SIMULINK MODEL 


%H%%%%%%%%%%%%UVH%MAMUVVVMMM%MMUMAVV%VVUV%VV%MAVVVAVHUVVVY%% 


% 

% 

% orbeom2d.m 

% 

% Equations of motion for 2d orbit 

% 

% Author: LT Jeffrey T. Elder 

% Naval Postgraduate School 
% November 1996 

% 


LUWPWLWLWLWLNWLWVVVA L/L /% A/V VULVA YAY NYY% 
% 


function xdot=orb_eom(x) 


mu = 1.327124399355el1; % Mass parameter for Sun 
AU = 1.49597870e8; % Astronomical unit 

r2 = x(3)*2 + x(4)2; % Radial distance squared 
mur2 = mu/r2; 

r =sqrt(r2); % Radial distance 

xdot(1) = -mur2*x(3)/r; % vx Ist ODE 

xdot(2) = -mur2*x(4)/r; % vy Ist ODE 

xdot(3) = x(1); % x Ist ODE 

xdot(4) = x(2); % y Ist ODE 
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APPENDIX D. ANALYTIC METHOD SCENARIO MODELS 


The following collection of figures represents the analysis performed using the 
numerical routine developed to perform the analytic solution. The scenarios are evaluated 
for impact true anomalies from -180° to +180° referenced to the asteroid orbit perihelion. 
The "steps" between impact scenarios is generally in 30° increments, however, for clarity 
purposes, the step size is reduced to 10° or even 1° intervals at times. The primary 30° 
increment figures occur where a "surface plot" and impact scenario plot appear together. 
At other increment values, two surface plots appear together. 

For the surface plots at impact scenarios approaching +180”, the surface appears 
rough due to the coarseness of the step size in impulse time and impulse direction. At 


smaller step sizes the behavior is smooth and well behaved. 
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